TNFSF15 and MIA Variant Associated with Immunotherapy and Prognostic Evaluation in Esophageal Cancer

Background Esophageal cancer (ESCA) is a common gastrointestinal tumor, and China is one of the regions with a high incidence. Tumor immune-related cells play important roles in the tumorigenesis and development of ESCA. However, the role of tumor immune-related genes in the development of ESCA has not been established. Methods In this study, weighted gene coexpression network analysis (WGCNA) was used to analyze ESCA gene expression using data from The Cancer Genome Atlas (TCGA) database. Gene expression was associated with clinical traits, and modules related to CD8+T cells, dendritic cells, and regulatory T cells (Tregs) were obtained. Results The GO analysis showed that inflammatory chemotaxis networks were activated by cell chemotaxis, chemokine activity, and chemokine binding receptor. Three hub genes (IL17C, TNFSF15, and MIA) related to tumor immunity and metastasis were identified by WGCNA, and the abnormal expression of each hub gene in ESCA has a poor prognosis, especially in patients with high expression (P < 0.05). The risk assessment analysis also showed that tumor stage was positively correlated with tumor risk in ESCA (P < 0.05). Therefore, more than 50 pairs of tumor tissues from the T1–T3 stages with different degrees of differentiation and paracancerous tissues were selected to confirm the expression of the three genes using RT-qPCR and immunofluorescence (IF). The infiltration of CD8+ T cells in tumor tissues was lower than that in normal tissues. According to the RT-qPCR, the expressions of IL17 C, TNFSF15, and MIA in moderately and poorly differentiated tissues were significantly higher than those in normal tissues (P < 0.05). In contrast, their expressions were decreased in high differentiated tissues (P < 0.05). Furthermore, IL17C, TNFSF15, and MIA were all positively correlated with immune checkpoint PD-1; TNFSF15 and MIA were also positively correlated with CTLA4, TIGIT, and CD96. Conclusion In summary, IL17C, TNFSF15, and MIA may act as biomarkers for prognosis in moderately and poorly differentiated ESCAs, and they may be used as predictive genes of immunotherapy associated with CD8+ T cell and Tregs invasion in ESCAs.


Background
Esophageal cancer (ESCA) is one of the most aggressive tumors and the sixth leading cause of cancer-related deaths globally. In developing countries, ESCA ranks eighth in incidence and ffth in mortality [1]. It is expected that the global morbidity and mortality of ESCA will increase within the next few decades [2]. Tere are two main pathological types of ESCA which are esophageal adenocarcinoma (EAC) and esophageal squamous cell carcinoma (ESCC). Simultaneously radiotherapy and chemotherapy are the main treatment options for advanced ESCA [3]. Due to the limitations of conventional therapies for ESCA, targeted immunotherapy is a key step in the individualized therapy of ESCA [4]. Immunotherapy targeted at programmed deathligand 1 (PD-L1) or programmed death (PD-1) checkpoints has been approved for ESCA [5,6]. However, events within the tumor immune microenvironment, including immune cell dysfunction or the accumulation of tumor metabolites, may afect the response of patients with ESCA to immune therapy [7]. As such, the prognosis is still very poor, especially for patients with advanced and metastatic ESCA. Terefore, exploring the mechanisms underlying ESCA immunity and identifying biomarkers is essential for early diagnosis, precision treatment, and prognostic evaluation of cancer.
Tumor-infltrating immunocytes have been implicated during diferent stages of tumorigenesis and development [8]. Activated CD8+ T cells play a pivotal role in tumor cytotoxicity, which releases some cytotoxic factors, leading to cell death [9]. Tregs and CD4 + FoxP3-T cells (TCONV) are capable of inhibiting antitumor immune responses, promoting tumorigenesis, and accelerating metastasis [10][11][12]. Tumor-infltrating immunocytes may afect tumor development and prognosis. In this study, we used WGCNA and GSEA to determine the prognostic markers associated with CD8+ T, CD4+ cells, and Tregs in ESCA. We identifed these molecules in 55 ESCA tissues and 55 paired normal adjacent tissues. Furthermore, the cBioPortal database, which includes the clinical and mutation data originating from 2201 patients, was used to analyze hub genes related to overall survival and mutation rates in ESCA. Tese results will be benefcial for our understanding of the role of CD4 + T cells in ESCA and provide a basis for ESCA treatment.

Gene Expression Data and Processing.
We obtained RNA expression and clinical data of ESCA from TCGA, including those for 165 tumor specimens and 11 normal specimens. First of all, we used the "org.Hs.eg.db" package in R language to convert the "Ensembl id" into a gene symbol (https:// www.bioconductor.org). At the same time, we converted RNA-seq-FPKM-count to RNA-seq-TPM-count. Finally, only 1078 immune-related genes downloaded from the ImmPort database (https://www.immport.org/) were selected for WGCNA.

Evaluation of Tumor-Infltrating Immune Cells (TIICs).
CIBERSORT, a general calculation method, was used to quantify the cellular part from the batch tissue gene expression profling (GEP). It can accurately estimate the immune components of the tumor biopsies [13]. Terefore, it was used to calculate the proportion of 22 kinds of TIICs in the study, and 10 types of TIICs that were diferent from normal tissues were selected. Te proportion of 5 T-cell subtypes in each sample was selected as the trait data for WGCNA.

GSEA Enriches the Immune-Related Pathways.
According to the median of gene expression, we divided the patients into two groups as follows: high and low expression groups. Te hallmark, c2, c5, and c7 gene sets downloaded from the Molecular Signatures Database (https://software. broadinstitute.org/gsea/msigdb/index.jsp) were enriched and analyzed using the GSEA_4.1.0. Te immune-related pathways were also enriched.

Construction of a Gene Coexpression Network for ESCA.
Te R software package "WGCNA" is used to construct a coexpression network with the expression values of 1078 genes to generate sample clusters and detect outliers [14]. "Flash-Clust" was used to screen the samples, and the characters of each sample were transformed into one color. Te criterion for selecting the soft threshold is that the degree of freedom is greater than 0.85; therefore, we set the soft threshold to four (scale-free R2 � 0.85) to construct the weighted adjacency matrix. Furthermore, the weighted adjacency matrix was transformed into a topological overlap matrix (TOM) to estimate the connectivity of genes in the network generation. Based on the diferential detection of TOM, we chose a cutting height of 0.25 and a minimum module size of 30 and divided the gene into nine diferent modules.

Determination of the Correlation between Modules and
Clinical Features. Module eigengenes (MEs) are defned as the principal components of each module. To determine the signifcance of the modules, we calculated the correlation between the MEs and the level of immune cell infltration in the tumor tissue using the Pearson test. At P < 0.05, an individual module was considered to be signifcantly related to immune cells. We selected the immune cell subtypes of interest and the module with the highest correlation coefcient and formed the hub module with them. In addition, we calculated the correlation between the genes and clinical features (cor. gene signifcance) and the correlation between the genes and MEs (cor. Gene module membership) to screen key genes in the hub module.

Multifactor Correlation Analysis and Pearson Correlation
Coefcient. We performed multivariate Cox proportional hazards regression analysis to determine the association of clinical parameters, including age, gender, tumor stage, and smoking, with a prognostic assessment of   ESCA. Risk scores were calculated using the following formula: In addition, the correlation between three genes and immune cells was, respectively, analyzed by the Pearson correlation coefcient. Furthermore, the same correlation coefcient was used to analyze the correlation between three hub genes and immune checkpoints. Pearson correlation coefcient was calculated using the following formula: (2)

Survival and Mutation Analyses.
Genes related to infammation and metabolism were selected as hub genes for the key genes. Next, we used Kaplan-Meier plots (https:// kmplot.com/analysis/) to analyze the correlation between the panoncogene (hub gene) mutation group and the normal control group in patients with ESCA. In addition, we obtained the mutation rate of hub genes from the cBioportal database (https://www.cbioportal.org/).

Patient
Samples. Fifty pairs of ESCC and adjacent normal tissues were obtained from the afliated Hospital of North Sichuan Medical College, including high diferentiation, middle diferentiation, and low diferentiation samples. All patients were surgical patients diagnosed with ESCA based on pathology. Te tumor stages were based on the eighth edition of the Union for International Cancer Control and American Joint Committee on Esophageal Carcinoma TNM classifcation system. Te clinical characteristics of these patients were obtained from their medical records. Te characteristics were age, sex, pathological tumor stage, and degree of tumor diferentiation degree (Table 1).

RNA Extraction and Real-Time Quantitative PCR (RT-qPCR).
Te tissues were lysed with Trizol reagent (TaKaRa), and the total RNA was extracted. Te concentrations of total RNA were detected using a NanoDrop 2000 C spectrophotometer (Termo Scientifc, USA). Te total RNA was reverse-transcribed into cDNA by Revertaid First Strand cDNA synthesis reagent (CAT: # K1622, Termo Fisher Science). Te cDNAs were amplifed by real-time quantitative PCR (RT-qPCR) with specifc primers (Table 2), and PowerUp SYBR Green Master Mix (CAT: # A25742, Termo Fisher Science) was used as the luminescent substrate. Finally, we analyzed the melting curve, and the expression of the hub gene was normalized with β-actin.

Infltration of Mature T and B Cells in ESCA Decreased
While Naïve B Cells Increased. Te analysis of adaptive immune cell infltration by 10 subtypes in 161 ESCA samples showed that the proportion of memory resting CD4 Tcells in the quiescent stage was the highest, followed by those of the CD8 T cells, memory activated CD4 T cells, follicular helper T cells, naïve B cells, and plasma cells (Figure 1(a)). Te diferences in immune infltrating cells between the normal and tumor samples were compared and analyzed, and we discovered that the degree of immune cell infltration differed between normal and tumor tissues. Te degree of infltration of CD8 T cells and plasma cells in the tumor was signifcantly lower than that in normal tissues (P < 0.05), while the infltration of naïve B-cells in the tumor was higher than that in the normal group (P < 0.05) (Figure 1(b)). As shown in Figure 1(c), the CD8 T cells were signifcantly reduced in tumors compared to paracancerous tissues.

Analysis of Immune-Related Pathways Using GSEA.
GSEA was used to analyze the diferences in tumor immune infltration and tumor metastasis-related pathways for the normal and tumor tissues. As shown in Figure 2(a), the degree of enrichment of characteristic genes of efector CD8+ T cells during the early stage of cancer was remarkably higher than that during the later stage (P < 0.05); the downregulated IL6 deprivation pathway genes were enriched in tumor tissues (P < 0.05) (Figure 2(b)), CD4 + FoxP3-T (TCONV) cell  characteristic genes that promote tumor migration and metastasis were enriched in tumors, while nature Treg cell characteristic genes were downregulated (Figure 2(c)), and melanoma metastasis characteristic genes promoting the tumor metastasis pathway were enriched in tumor tissues (Figure 2(d)).

Weighted Coexpression Network of ESCA.
We obtained 1100 gene expression matrices and clinical information related to adaptive immunity in ESCA tissues from the TCGA database, excluding one outlier sample, and constructed a hierarchical clustering tree (Figure 3(a)) for the remaining 76 samples. In addition, we labeled the immune cells with the greatest infltration in cancer tissues according to the immune landscape in the ESCA. Besides, we observed signifcant diferences from normal tissues under the resulting tree (Figure 3(b)). Following the proposal of the fle package, the soft threshold power value was used to construct a gene coexpression network (Figure 3(c)). Based on the construction of the coexpression module clustering and hierarchical clustering trees, a total of nine modules were identifed (Figure 3(d)).

Correlation between Modules and Adaptive Immune Cells with Signifcant Diferences in Tumor Tissues.
We determined the coexpression genes of each module according to their gene characteristics and calculated the correlation of the characteristic value of each module based on the characteristics of the adaptive immune cells, the cells with higher correlation are CD8 + T cells, naïve B cells, plasma cells, memory resting CD4 T cells, and regulatory T cells. We selected two modules with the highest correlation in the CD8 + T cells, as denoted by the ME-turquoise (R2 � 0.38, P < 0.05) and ME-blue rows (R2 � 0.29, P < 0.05). At the same time, we selected two modules with the highest correlation in the regulatory T cells, as denoted by the MEturquoise (R2 � 0.33, P < 0.05) and ME-brown (R2 � 0.59, P < 0.05) rows (Figure 4). Te bubble chart shows similar results (Supplementary Figure 1).

Tree Hub Genes Were Genetic Variants of ESCA.
Based on WGCNA analysis, the expression of three hub genes IL17C, TNFSF15, and MIA were all negatively associated with dendritic cells. Among them, TNFSF15 expression was also negatively associated with CD4 + memoryactivated T cells. In addition, TNFSF15 and MIA expression were positively associated with Tregs and CD4 + memory resting T cells, respectively (P < 0.05) ( Figure 5(a)). In addition, the genetic variability of these genes in ESCA was analyzed. Te analysis of genetic variation in the cBioPortal database showed that the proportion of IL17C, TNFSF15, and MIA in ESCA was 2.0%, 1.7%, and 1.9%, respectively (P < 0.05), and they all presented several alterations including amplifcation and deletion in patients with ESCA ( Figure 5(b)).

TNFSF15, IL17C, and MIA Are Important Biomarkers for
Evaluating the Prognosis of Patients with ESCA. We further evaluated the value of the above hub gene selection in evaluating the prognosis of patients with ESCA, and the overall survival (OS) of the patients with the three abnormal hub genes was analyzed using the cBioPortal database. Te cBioPortal database contains 2201 clinical samples from eight ESCA-related research websites. As shown in Figure 6(a)-6(c), the patients with ESCA with abnormalities of IL17C, TNFSF15, and MIA had remarkably lower 5-year and 10-year survival rates than the controls, and their overall survival duration was not more than fve years (P < 0.05). Te overall survival duration of patients with IL17C gene abnormalities was shorter, even less than one year (P < 0.05). Moreover, we also analyzed the relationship between the expression of the three genes and prognostic survival using the TCGA database. Te results showed that patients with low expression of IL17C, TNFSF15, or MIA lived longer than those patients with high expression of one of them (P < 0.05) (Figure 6(d)-6(f )).

IL17C, TNFSF15, and MIA Expression Correlated with the Degree of Diferentiation of ESCA, Validated in Tissues from Stage T1 to T3 Esophageal Squamous Cell Carcinoma.
Some clinical parameters: age, gender, stage, and smoking are prognosis evaluations in the clinical. As shown in Figure 7(a), there was a positive correlation between the ESCA stage and prognosis evaluation (P < 0.05). Terefore, tissues from patients with stage-T1 to T4 ESCA who were only treated with surgery were selected to verify these candidate genes. Te tumor stage was determined based on the results of pathological HE staining (Figures 7(b) and 7(c)). Ten, the qPCR assay was performed for 55 pairs of clinical samples, including adjacent and tumor tissues. As shown in Figure 7, the expressions of IL17C, TNFSF15, and MIA in tumor tissues were signifcantly diferent from those in adjacent tissues. In tissues with moderate and low diferentiation, the expressions of these three genes were signifcantly higher than those in tumor-adjacent tissues (P < 0.05) (Figures 7(d)-7(f )). In highly diferentiated ESCA tissues, by contrast, the expressions of IL17C, TNFSF15, and MIA  Journal of Oncology in tumor tissues were signifcantly lower than those in normal tissues (P < 0.05) (Figures 7(g)-7(i)).

Relationship between the Tree Genes and Immune
Checkpoints in ESCA. PD-1, CTLA4, CD96, and TIGIT are important immune checkpoints that are potential targets for tumor immunotherapy. To further understand the role of these genes in ESCA, the relationship of the three genes with those immune checkpoints was assessed. As suggested in Figure 8, MIA and TNFSF15 expression were signifcantly positively correlated with PD-1, CTLA4, CD96, and TIGIT in ESCA (P < 0.05); IL17C expression was positively correlated with PD-1 (P < 0.05) (Figure 8).

Discussion
CD8+ T cells play an essential role in efective immune responses. In our study, the infltration of cytotoxic T cells was decreased in tumor tissues, suggesting that tumor immunity was disregulated. Furthermore, CD4 + FoxP3-T (TCONV) cell characteristic genes that promote tumor migration and metastasis were enriched in tumors, while nature Treg cell characteristic genes were downregulated. Terefore, to determine the characteristics of adaptive immunity-related genes in ESCA, we performed RNA-seq analysis of the ESCA samples and compared the expressions of adaptive immunityrelated genes with that of normal tissues, and three modules related to CD8 T cells and Tregs were obtained. On this basis, IL17C, TNFSF15, and MIA which were related to tumor immunity, invasion, and metastasis were identifed. IL17C, an autocrine cytokine, is an important factor in the innate immunity of the epithelium. IL17C may promote or inhibit tumorigenesis by regulating immune cell function [15][16][17]. In our clinical samples, it was amplifed in moderately and poorly diferentiated tissues but downregulated in high diferentiated tissues. Te survival duration of patients with altered IL17C was less than 30 months, and patients with high expression of that were living shorter than those with low expression. TNFSF15 belongs to the tumor necrosis factor (TNF) ligand family induced by TNF and IL-1α, which plays an important role under disease conditions, such as cancer and stroke, by maintaining vascular and lymphatic vessel homeostasis [18][19][20][21]. In our study, the expression of TNFSF15 was signifcantly higher in moderately and poorly diferentiated tissues than in normal tissues, and patients with TNFSF15 alterations were in poor prognosis. All this suggested that TNFSF15 may play a risk role in moderate and low diferential ESCA. MIA is a neurotransmitter involved in tumor invasion and migration. MIA is a melanomaderived growth regulatory protein that can inhibit the growth of melanoma cells and some other neuroectodermal tumors (including gliomas) in vitro, and it may be an antitumor molecule in melanoma [22]. But in ESCAs, the expression of MIA was increased in moderately or poorly diferentiated tumor tissues, and patients with high diferentiated tumors lived longer than that middle-low diferentiated patients, which indicates that MIA may have diferent roles according to the diferent tumor. All this suggests that IL17C, TNFSF15, and MIA may be prognostic risk factors for OS in middle-low diferentiated ESCA, respectively. Furthermore, the analysis of immunocyte infltration in ESCA suggested that CD8+ T cells were signifcantly lower than that in paracancerous tissue, and we used IF to confrm that CD8+ T cells were reduced in ESCA tissues. PD-1, CTLA4, CD96, and TIGIT are associated with each other in tumor immunity, which are candidates for immunotherapy [23][24][25][26]. IL17C and TNFSF15 are involved in tumorigenesis and development by impacting immune cell function, and MIA is associated with tumor invasion [27][28][29]. Ten, we analyzed the relations among IL17C, TNFSF15, and MIA and the immune checkpoint. TNFSF15 and MIA were positively associated with the four immune check-points, PD-1, CTLA4, CD96, and TIGIT. In addition, IL17C was only positively associated with PD-1 expression, and we used IF to confrm that IL17C was highly expressed in ESCA tissue. Tese results indicate that highly expressed IL17C, TNFSF15, and MIA in ESCA may infuence the development and prognosis of ESCA through participating in the immune regulation process, which may act as predictors of immunotherapy for ESCA.
However, there were some limitations in our study, such as the sample size limiting the accuracy of the analysis and the experimental results. Furthermore, more time was needed to review patients with special gene phenotypes and evaluate their prognosis. Finally, the specifc mechanisms underlying the actions of these alteration genes in ESCA need to be elucidated in future research.
In summary, we used multiple biological information analyses and experiments involving clinical samples to identify potential biomarkers of ESCA. Tree hub genes were identifed to be abnormally expressed in tumors, and they afected the prognosis of patients with ESCA. Te immune cell function regulating IL17C, TNFSF15, and MIA have been identifed as potential biomarkers for the evaluation of ESCA diagnosis and prognosis, which may also be predictors of immunotherapy.

Conclusions
In this study, CD8+ T cells were decreased in tumor tissues which contributed to the immune imbalance of the tumor. WGCNA was used to identify the characteristics of immunerelated genes in ESCA. Tree modules related to CD8+ T cells and Tregs were identifed, and three hub genes related to the survival duration and prognosis of ESCA were obtained, which are IL17C, TNFSF15, and MIA, and they can be used as therapeutic targets for further study and clinical markers to predict the prognosis of patients.

Data Availability
Publicly available datasets were used in this study. Te data can be found in the TCGA database and the cBioportal website.

Ethical Approval
Tis study was approved by the Ethics Committee of the Afliated Hospital of North Sichuan Medical College.

Consent
Informed consent was obtained from each patient at the beginning of the study.

Conflicts of Interest
Te authors declare no conficts of interest.